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ABSTRACT 

This  work  has  involved  the  use  of  the  Naval  Postgraduate  School's  TRITECH  ST  725 
high  frequency  mechanically  scanned  sonar  system  to  acquire  sonar  images  of  simulated 
surface  roughness  on  an  aluminum  plate.  Signal  post  processing  for  such  image  data  is 
reviewed,  and  processed  data  is  analyzed  and  compared  to  the  known  roughness  locations 
on  the  plate. 

The  simulated  roughness  (a  pattern  of  one  half  inch  steel  nuts)  is  used  in  a  preliminary 
experiment  as  part  of  the  development  of  a  sonar  detection  system  for  marine  growth  on 
ship  hull  plating.  Such  a  sonar  system  will  be  an  integral  part  of  any  Ship  Hull 
Autonomous  Cleaning  Robot  (SHACR). 

Contained  in  this  report  is  a  description  of  the  experimental  arrangement,  typical  sonar 
returns,  a  summary  of  image  processing  techniques  appropriate  to  this  problem,  and  results 
of  processed  data  as  compared  to  the  known  locations  of  the  simulated  roughness. 

The  algorithms  presented  here  will  ultimately  lead  to  a  real  time  processing  capability  for 
the  specification  of  location,  extent,  roughness  level,  as  needed  for  the  automatic  direction 
of  a  SHACR's  motion  control  and  cleaning  systems. 
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1.0    INTRODUCTION 

Ship  hull  fouling  is  a  common  phenomenon  having  a  strong  impact  on  the  Navy's 
operating  costs.  Fouling  takes  place  through  the  growth  of  marine  life,  as  in  barnacles  and 
grass,  that  increase  the  drag  on  the  ship  while  it  is  motion.  There  is  an  attendant  loss  of 
efficiency  and  an  increase  in  fuel  consumption  costs.  On  the  other  hand,  cleaning  fouled 
hulls  is  also  expensive,  and  may  require  dry  docking.  In-situ  cleaning  without  dry  docking 
carries  the  risk  of  cleaning  too  deeply  and  releasing  environmentally  damaging  deposits 
(and  particularly  copper)  into  the  harbor. 

A  potentially  cost  saving,  environmentally  friendly,  method  of  cleaning  ship  hulls 
involves  the  use  of  underwater  robots  equipped  with  sensors  to  detect  and  locate  fouled 
sections  of  ship  hull  plating  thereby  selecting  only  the  appropriate  tools  for  the  cleaning  job 
which  could  change  from  area  to  area.  Such  a  device  would  be  useful  as  an  aid  to  divers 
who  currently  have  a  difficult  time  in  clearly  identifying  the  location  and  extent  of  marine 
fouling  underwater  -  even  with  the  use  of  video  cameras. 

The  purpose  of  this  study  was  to  initiate  the  development  of  sonar  based  sensory 
techniques  for  use  with  an  autonomous  underwater  vehicle  that  will  be 

(1)  able  to  better  detect  the  extent  of  fouling  on  the  hull  of  the  ship,  and, 

(2)  localize  the  areas  of  significant  fouling 

allowing  selective  cleaning  to  be  put  into  effect. 

In  discussions  with  the  sponsor,  it  was  established  early  that  visual  sensing  through 
video  cameras  would  not  be  a  generally  useful  technology  because  the  murky  waters  of  a 
harbor  would  prevent  video  images  from  being  clearer  than  the  current  capability  offers 
using  divers. 

The  study  has  therefore  focused  on  the  use  of  high  frequency  imaging  sonar  as  a 
device  to  detect  and  localize  segments  of  fouling.  In  the  following  sections  of  this  report, 
we  will  describe  the  particular  sonar  used,  the  example  of  an  aluminum  plate  with  simulated 
roughness,  and  the  signal  processing  techniques  available  for  use  in  converting  data 
streams  from  the  sonar  head  into  signals  that  could  be  used  to  direct  the  movement  of  the 
cleaning  robot. 

While  the  study  was  funded  at  a  very  low  level,  the  results  do  indicate  that  further 
investigations  are  warranted  with  the  potential  of  reasonable  success  on  real  hulls.  We 
strongly  recommend  that  further  experiments  be  conducted  on  plates  that  have  realistic 
examples  of  marine  fouling,  followed  by  an  experimental  program  to  investigate  hulls  in 
situ. 
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Figure  1  Sketch  of  the  NPS  AUV  II  Autonomous  Underwater  Vehicle 


2.0   THE  TRITECH  ST  725  SONAR 

2. 1       General 

The  use  of  high  frequency  sonars  (to  1  Megahertz)  is  expected  to  provide  usable 
results  for  the  imaging  of  marine  fouling  from  barnacles  and  grass. 

At  the  Naval  Postgraduate  School  we  have  two  TRITECH  high  frequency  sonars. 
These  are  small  units  that  fit  into  a  small  autonomous  underwater  vehicle  -  the  NPS  AUV  II 
-  shown  in  Figure  1.  They  would  be  suitable  for  installation  on  a  SHACR,  and  this  is  one 
reason  why  interest  in  their  use  is  strong. 

The  TRITECH  ST  725  is  a  sector  scanning  sonar  that  is  mechanically  scannable 
through  full  360  degrees  rotation  -  or  may  be  swept  through  any  arbitrary  defined  sector 
given  by  starting  and  terminating  angles.  In  sweeping  through  a  defined  sector,  the  sonar 
steps  at  0.9  degrees  per  step  at  a  rate  of  20  steps  per  second,  where  at  each  step  (bearing),  a 
ping  is  emitted  and  the  return  is  converted  into  an  intensity  (average  power)  value  contained 
in  each  of  several  discrete  segments  of  range  called  range  bins.  There  are  normally  128 
range  bins  associated  with  the  sonar  returns  from  the  ST  725  unit.  The  effective  range 
associated  with  any  bin  is  therefore  1/128  of  the  maximum  range  setting  in  use.  One 
exception  to  this  is  when  the  maximum  range  setting  of  6  m  is  used  where  it  is  divided  into 
only  64  bins  associated  with  600/64,  or  a  discretization  of  9.375  cm.  for  the  shortest  range 
setting.  It  follows  that  for  each  bearing  setting,  a  "scan  line"  is  developed  that  provides  a 
vector  of  intensity  values  -  one  for  each  range  bin.  The  image  processing  problem  is  to 
convert  the  data  so  generated  into  meaningful  metrics  as  the  sonar  sweeps  over  the  plating 
of  the  ship  hull.  The  metrics  defined  must  be  such  that  they  will  indicate  the  location  and 
general  degree  of  fouling. 

The  second  sonar,  the  TRITECH  ST  10(H)  is  a  one  megahertz  sonar  with  a  one 
degree  conical  beam,  mechanically  scanned,  where  range  for  each  bearing  is  based  on  time 
of  flight  of  the  first  strong  return.  The  data  received  at  each  bearing  is  therefore  directly 
related  to  the  distance  from  the  sonar  head  to  the  object  of  the  first  return,  and  the  ST  1000 
thus  becomes  a  useful  profiling  device. 

As  an  example,  Figure  2  shows  a  typical  result  of  processing  returns  from  the  ST 
725  sonar  being  set  for  a  maximum  range  of  6  meters  in  the  NPS  Hovering  Tank  -  a  twenty 
foot  square  tank  where  the  tank  walls  are  clearly  identified.  In  the  basic  configuration,  the 
ST  725  communicates  with  a  surface  located  PC  computer  with  a  VGA  graphics  screen  to 
display  the  intensity  results  on  a  polar  plot.  The  intensity  of  each  pixel  (  a  pixel  is  defined 
as  the  area  of  the  screen  occupied  by  a  single  range  bin  at  a  particular  bearing  step),  is 
plotted  in  a  color  selected  according  to  a  color  /  intensity  map. 

However,  while  it  is  easy  to  see  the  sonar  image  on  the  PC  screen,  this  mode  of 
processing  is  unsuitable  for  the  integration  into  a  SHACR  control  system  architecture. 
What  is  needed  is  the  automatic  identification  of  location  and  degree  of  fouled  segments  of 
the  hull  plating.  The  algorithms  developed  for  that  purpose  would  be  the  'methods'  of  a 
Tactical  Level  Object  in  the  context  of  the  Tri-Level  Software  Architecture  of  the  RBM,  -  a 
mission  package  controller  -  as  proposed  for  the  NPS  AUV  II,  but  which  could  be  easily 
adapted  to  a  SHACR.  A  large  amount  of  work  is  still  necessary  to  define  the  best 
methods  which  may  be  used  in  raising  the  level  of  autonomy  in  hull  cleaning  robots, 
although  the  results  of  studies  of  the  kind  performed  here  will  certainly  help. 
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Figure  2  Typical  Sonar  Image  Display  from  the  ST  725  in  the  NPS  Hover  Tank 
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2 . 2  Modes  of  Operation 

The  ST  725  sonar  can  be  set  to  have  a  maximum  range  at  either  6,  10,  20,  25,  50, 
75,  or  100  meters.  A  gain  setting  is  adjustable  by  the  operator,  and  can  be  set  so  high  that 
all  returns  are  of  the  maximum  intensity  producing  meaningless  results  if  care  in  its  setting 
has  not  been  taken.  For  an  appropriate  gain  setting  -  depending  on  the  distance  and 
conditions  of  the  targeting  -  distinction  can  be  made  between  returns  from  roughness  and 
smooth  segments  of  plating.  The  intensity  of  the  returns  are  digitized  into  a  scale  of  1  -  15, 
where  15  represents  the  highest  intensity  of  return.  The  sonar  displays  can  be  thresholded 
and  controlled  manually  prior  to  display  and  the  sweep  speed  can  be  either  at  0.9  degrees 
per  step,  or,  for  low  resolution  but  faster  analysis,  at  1.8  degrees  per  step. 

2.3  Automatic  Control  of  the  ST  725  Sonar 

We  have  been  working  -  separately  from  this  project  -  to  automatically  control  the 
gain,  threshold,  sweeping  sector,  and  sweep  rate  through  ASCII  commands  sent  to  the 
serial  port  of  the  ST  725  sonar.  Automatic  control  of  the  initialization  process  for  the  sonar 
has  also  been  completed  so  that  separate  code  -  part  of  the  execution  level  of  the  NPS  AUV 
II  control  code  will  be  able  to  direct  the  settings  of  the  sonar  and  be  able  to  read  and  process 
the  returns.  This  work  is  soon  to  be  available  as  a  useful  real  time  controller  for  the  ST 
series  of  sonars. 


3.0    IMAGE  ACQUISITION  CONCEPTS 


3. 1  General 

Not  the  least  part  of  the  image  acquisition  problem  is  the  arrangement  of  the  sonar 
head  in  relation  to  the  plate  whose  roughness  characteristics  are  to  be  measured.  Potentially 
there  are  a  large  number  of  configurations  that  could  be  used  to  acquire  data  representative 
of  fouling.  Two  concepts  have  been  considered. 

For  instance,  a  surface  roughness  profile  could  be  measured  using  the  ST  1000  as  a 
vertical  profilometer  where  the  roughness  of  the  plate  surface  would  be  the  measure  of  the 
degree  of  fouling.  Because  of  the  resolution  obtained  with  the  least  significant  bit  at  the 
smallest  range  setting,  the  ST  1000  sonar  would  be  expected  to  be  accurate  to  about  2  cm. 
of  height  change  .  This  was  considered  too  coarse  to  be  of  primary  interest  to  this 
application  and  it  is  concluded  that  a  direct  surface  roughness  measurement  using  this  sonar 
is  not  sufficiently  precise  unless  a  special  purpose  sonar  is  developed. 

Alternatively,  using  the  ST  725  as  a  sector  scanner  the  24  degree  high  by  1  degree 
wide  beam  could  be  used  to  ensonify  an  area  of  the  plate  as  it  sweeps  horizontally  across 
the  surface.  Returns  from  the  ensonified  regions  would  then  be  distinguished  by  their 
intensities.  While  the  resolution  from  the  ST  725  would  create  pixels  a  little  larger  than  the 
ultimately  desired  size  (an  unknown  quantity  at  the  present  time)  the  horizontal  arrangement 
appears  to  offer  the  most  promise.  The  geometry  for  such  a  setup  is  illustrated  in  Figure  3. 

3.2  Horizontal  Configuration  Detail 

The  side  view  in  Figure  3  shows  the  sonar  beam  ensonifying  a  section  of  the  plate 
with  some  roughness  as  modeled  by  the  circles  together  with  a  return  ray  from  the  second 
object  in  the  range  bin  marked  as  No.  2. 

At  the  nominal  range  setting  of  6  meters,  there  are  64  range  bins  and  so  the  length 
of  each  range  bin  is 

6/128  =  9.375  cm.  per  bin 

For  a  beam  width  of  1  degree,  the  beam  width  becomes 


Bin  Number 

Range  cm. 

Width  at  End  of  Bin  cm. 

1 

00.000  -  09.375 

0.1636 

2 

09.375  -  18.750 

0.3272 

3 

18.750-28.125 

0.4909 

4 

28.125 -37.5(H) 

0.6545 

5 

37.500  -  46.875 

0.8181 

6 

46.875  -  56.25 

0.9817 

7 

56.250  -  65.625 

1.1453 

8 

65.625  -  75.000 

1.3089 

9 

75.000  -  93.750 

1.4726 

Each  sonar  pixel  is  thus  of  radial  length  equal  to  9.375  cm  with  a  tangential  width 
equal  to  the  arc  subtended  by  one  degree  at  the  appropriate  radius  from  the  sonar  head. 
Average  intensity  of  return  from  each  sonar  pixel  is  what  is  measured  by  the  sonar 

3.3  ST  725  Default  Data  Acquisition 

The  ST  725  Sonar  is  equipped  with  a  software  executable  file  running  on  an  MS- 
DOS  PC  with  a  VGA  Graphics  screen  capability  that  will  process  the  returns  from  the 
sonar  in  real  time  -  scanline  by  scanline  and  present  the  plots  of  intensity  on  a  polar  plot 
indicating  Cartesian  X- Y  coordinates  of  the  returns  plotted.  While  this  mode  of  operation 
is  satisfactory  for  human  operator  visualization  of  the  sonar  returns,  it  is  not  suitable  for 
consumption  in  real  time  by  a  SHACR. 

3.4  Typical  Data  File  Format 

The  sonar  responds  to  commands  that  are  issued  through  its  serial  port.  A  listing  of 
possible  commands  includes: 

1 .  Ping  and  process  returns 

2.  Rotate  one  step  (clockwise  or  counterclockwise) 

3 .  Read  serial  port  for  processed  return  data 

4.  Initialize  head  to  default  settings 

5.  Set  range,  gain,  resolution  data 

To  generate  an  image,  the  sonar  code  that  operates  the  head  issues  a  series  of 
commands  that  request  the  sonar  to  ping,  process  the  returns,  and  take  the  next  step.  At 
this  point,  the  returns  data  are  written  to  a  file  which  the  user  id  asked  to  name,  and 
displayed  on  the  graphics  screen  for  visualization. 

Examination  of  the  data  file  format  has  revealed  that  it  is  typically  a  large  file  in 
which  -  with  some  degree  of  compression  by  eliminating  all  data  with  intensity  values 
equal  to  zero  -  the  format  of  the  data  is  as  shown  in  the  Table  1  below.  What  is  clear  from 
the  form  of  the  data  in  the  Table  1  above  is  that  each  scanline  is  associated  with  a  particular 
bearing  angle.  Within  each  scanline,  non  zero  intensity  readings  are  associated  with  range 
bins  at  various  range  values.  By  reading  the  PC  file  created  by  the  default  program,  we  are 
able  to  store  such  data  and  replot  using  the  VAX  system  at  NPS  to  produce  a  replica  of  the 
screen  data  produced  by  the  default  data  processing.  Replays  of  data  taken  in  the  NPS 
hovering  tank  where  the  ST  725  ensonifies  the  aluminum  plate  described  later,  are  now 
given. 

Figure  4  shows  a  typical  sonar  intensity  map  corresponding  to  a  single  sweep  over 
a  rough  aluminum  plate.  The  data  is  taken  from  a  file  created  with  a  sonar  gain  setting  of 
1 1  and  the  data  is  plotted  covering  2000  lines  of  data  giving  a  sector  sweep  of  an  included 
angle  of  90  degrees. 

High  points  of  intensity  in  the  returns  are  shown  in  black  while  lower  intensity 
returns  are  shown  in  colors  nearer  to  white  which  corresponds  to  zero  intensity.  The  series 
of  figures  in  Figure  4  show  the  effects  of  gain  setting  and  later,  the  effect  of  a  low  versus 
high  resolution. 


TABLE  1 


Table  of  Typical  Data 
Data  Line  No.  Bearing  Degrees  Range  (cm) 


-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
•39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
-39.600000 
2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 

2  -38.700000 


Intensity  Measure 


0.000000 

15 

0.093750 

9 

0.187500 

3 

0.468750 

8 

0.562500 

1 

0.937500 

1 

1.031250 

5 

1.218750 

8 

1.312500 

3 

1.500000 

6 

1.593750 

4 

1.687500 

7 

1.781250 

8 

1.875000 

6 

1.968750 

10 

2.062500 

2 

2.156250 

1 

2.250000 

3 

2.343750 

5 

2.812500 

3 

2.906250 

1 

3.750000 

6 

3.843750 

7 

4.218750 

2 

4.312500 

10 

4.406250 

2 

4.781250 

4 

4.968750 

8 

5.062500 

4 

5.250000 

9 

5.343750 

4 

5.812500 

2 

0.000000 

15 

0.093750 

7 

0.187500 

2 

0.468750 

3 

0.937500 

1 

1.031250 

6 

1.125000 

5 

1.218750 

7 

1.312500 

6 

1.406250 

8 

1.500000 
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Figure  4c  Typical  Sonar  Intensity  Map 
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4.  SONAR  SIGNAL  PROCESSING 

4.1  Introduction 

Sonar  signal  processing  may  be  subdivided  into  a  four  principal  areas:  (a)  Preprocess- 
ing, (b)  Segmentation,  (c)  Description,  and  (d)  Interpretation.  Preprocessing  deals  with 
techniques  such  as  noise  reduction  and  enhancement  of  details.  Segmentation  is  the  process 
that  partitions  an  image  into  objects  of  interest.  Description  deals  with  the  computation  of 
features  like  size,  shape,  etc.,  suitable  for  differentiating  one  object  from  another.  Finally, 
interpretation  assigns  meaning  to  an  ensemble  of  recognized  objects.  In  our  case  where  we 
take  sonar  images  of  the  fouled  hull  of  a  ship,  the  data  will  be  interpreted  to  assess  the 
location  and  extent  of  fouling.  In  the  next  sections  we  describe  in  brief  the  methods  of 
preprocessing,  segmentation,  description,  and  interpretation  of  sonar  images. 

4.2  Preprocessing  of  data  from  the  sonar 

Preprocessing  may  be  considered  as  a  low-level  imaging  process  that  is  primitive  in  the 
sense  that  it  requires  no  intelligence.  It  will  take  us  from  the  image  formation  process  to 
the  compensations  such  as  noise  reduction,  and  finally  to  the  extraction  of  primitive  image 
features  such  as  intensity  discontinuities. 

There  are  two  basic  approaches  to  preprocessing:  the  spatial  domain  methods,  and  the 
frequency  domain  methods  via  the  Fourier  transform.  The  spatial  domain  refers  to  the 
aggregate  of  pixels  composing  an  image  and  the  spatial  domain  methods  are  procedures 
that  operate  directly  on  these  pixels.  Preprocessing  functions  in  the  spatial  domain  may  be 
expressed  as 

g(x,y)  =  h[f(x,y)}  (1) 

where,  f{x,y)  is  the  input  image,  g{x,  y)  is  the  preprocessed  image,  and  h  is  an  operator 
on  /,  defined  over  some  neighborhood  of  (x,y).  The  principal  approach  used  in  defining  a 
neighborhood  about  (x,y)  is  to  use  a  square  or  rectangular  subimage  area  centered  at  [x,y) 
as  shown  in  Fig.4. 

One  of  the  spatial  domain  methods  used  most  frequently  is  based  on  the  use  of  so-called 
convolution  masks  (also  referred  to  as  templates,  windows,  or  filters).  Basically,  a  mask  is 
a  2-D  array  (e.g.  3x3),  such  as  the  one  shown  in  Fig.4,  whose  coefficients  are  chosen  to 
detect  a  given  property  in  an  image.  As  shown  in  Fig. 5,  if  we  let  w\ ,  w2,  ■  ■  ■ ,  wq  represent  mask 
coefficients  and  consider  the  8  neighbours  of  (x,  y),  we  may  perform  the  following  operation 
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h  [f(x,  y)]  =w1f(x-l,y-\)  +  w2  f{x  -  l,y)  +  w3  f(x  -  1,  y  +  1) 
+w4f{x,y-  1)  +  w5  f(x,y)  +  w6  f(x,y  +  1) 
+w7f(x+l,y-l)  +  wsf(x+l,y)  +  wgf(x  +  l,y  +  l)  (2) 

The  frequency  domain  refers  to  the  aggregate  of  complex  pixels  resulting  from  taking 
the  Fourier  transform  of  an  image.  Due  to  the  extensive  processing  requirements,  frequency 
domain  methods  are  not  nearly  as  widely  used  as  are  spatial  domain  methods. 

In  the  next  subsections  we  discuss  the  methods  of  image  preprocessing  that  we  will  use 
in  the  next  part  of  this  project. 

4.2.1  Smoothing 

Smoothing  operations  are  used  for  reducing  noise  and  other  spurious  effects  that  may 
be  present  in  an  image  as  a  result  of  sampling,  quantization,  transmission,  or  disturbances 
in  the  environment  during  image  acquisition.  In  our  case,  we  have  observed  that  the  data 
acquired  by  the  sonars  have  a  significant  amount  of  noise  that  may  have  resulted  from 
multiple  reflections  of  the  sonar  signal. 

Neighborhood  Averaging  Neighborhood  averaging  is  a  straightforward  spatial  do- 
main technique  used  for  image  smoothing.  Given  an  image  f(x,y).  the  procedure  is  to 
generate  a  smooth  image  g{x,y)  whose  intensity  at  every  point  (x,y)  is  obtained  by  averag- 
ing the  intensity  values  of  the  pixels  of  /  contained  in  a  predefined  neighborhood  of  {x,y). 
In  other  words,  the  smoothed  image  is  obtained  by  using  the  relation 

<Ax,y)  =  ±     Y,      ^n'm)  ^ 

(n,m)6S 

for  all  x  and  y  in  f(x,y).  S  is  the  set  of  coordinates  of  all  points  in  the  neighborhood  of 
(x,  y),  including  (x,y)  itself,  and  P  is  the  total  number  of  points  in  the  neighborhood. 

Median  Filtering  One  of  the  principal  difficulties  in  neighborhood  averaging  is  that  it 
blurs  images  and  other  sharp  details.  The  blurring  can  often  be  reduced  significantly  by  the 
use  of  so-called  median  filters,  in  which  we  replace  the  intensity  of  each  pixel  by  the  median 
of  the  intensities  in  a  predefined  neighborhood  of  that  pixel,  instead  of  by  the  average.  A 
principal  function  of  the  median  filtering  is  to  force  functions  with  very  distinct  intensities 
to  be  more  like  their  neighbors,  thus  actually  eliminating  spikes  that  appear  isolated  in  the 
area  of  the  filter  mask. 

Image  Averaging  Consider  a  noisy  image  g{x,y)  which  is  formed  by  the  addition  of 
noise  n{x,y)  to  an  uncorrupted  image  f(x,y),  that  is 
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9{x,y)  =  f(x,y)  +  n(x,y)  (4) 

where,  it  is  assumed  that  the  noise  is  uncorrelated  and  has  zero  average  value.  The  objective 
of  this  procedure  is  to  obtain  a  smoothed  result  by  adding  a  given  set  of  noisy  images, 
gt(x,y),i  =  1,2,  •■•, K.  If  the  noise  satisfies  the  assumptions  stated  above,  the  it  can  be 
shown  that  if  the  image  g(x,y)  is  formed  by  averaging  K  different  noisy  images, 

1    K 
9(x,y)  = —^gi(x,y)  (5) 

t=i 
it  follows  that 

E{g(x,y)}  =  f(x,y)  (6) 

and 

°j{x,y)  =  J£  <*l(x,v)  (") 

where,  E{g{x,y)}  is  the  expected  value  of  g,  and  oUx,y)  and  cr*(x,y)  are  the  variances  of 
g  and  n,  all  at  coordinates  (x,y).  Since  E{y{x,y)}  =  f(x,y),  this  means  that  g(x,y)  will 
approach  the  uncorrupted  image  f{x,y)  as  the  number  of  noisy  images  used  in  the  averaging 
process  increases. 

4.2.2  Enhancement 

Image  enhancement  is  a  major  area  in  image  processing  and  much  of  the  success  of  the 
subsequent  processing  algorithms  is  dependent  upon  it.  Our  ability  to  accurately  predict 
the  areas  of  fouling  on  the  hull  of  the  ship  and  the  extent  of  the  fouling  will  depend  upon 
the  enhancement  procedures. 

Histogram  Equalization  Let  the  variable  r  represent  the  intensity  of  pixels  in  an  image 
to  be  enhanced.  It  will  be  assumed  initially  that  r  is  a  normalized,  discrete  variable  lying  in 
the  range  0  <  r  <  1.  For  r  lying  in  this  range,  attention  will  be  focussed  on  transformations 
of  the  form 

s  =  T(r)  (8) 

which  produce  an  intensity  value  s  for  every  pixel  value  r  in  the  input  image.  It  is  assumed 
that  the  transformation  satisfies  the  conditions: 

1.  T(r)  is  single- valued  and  monotonically  increasing  in  the  interval  0  <  T{r)  <  l. 

2.  0  <  T{r)  <  1  for  0<r<  1. 


Condition  1  preserves  the  order  from  black  to  white  in  the  intensity  scale,  and  condition  2 
guarantees  a  mapping  that  is  consistent  with  the  allowed  0  to  1  range  of  pixel  values.  The 
probabilities  are  given  by  the  relations 

Pr(rk)  =  — ,  0<r<l 

n 

fc  =  0,l,--,L-l  (9) 

where  L  is  thenumber  of  discrete  intensity  levels,  Mrfc)  is  an  estimate  of  the  probability  of 
intensity  rk,  nk  is  the  number  of  times  this  intensity  appears  in  the  image,  and  n  is  the  total 
number  of  pixels  in  the  image.  A  plot  of  pr(rk)  versus  rk  is  usually  called  a  histogram,  and 
the  technique  used  for  obtaining  a  uniform  histogram  is  known  as  histogram  equalization  of 
histogram  linearization.  Equation  (8)  in  the  discrete  form  is  given  by 


<      '     n 


n 

k 

=  J>-fo)  (10) 

3  =  0 

for  0  <  rk  <  1.  and  k  =  0, 1,  ■  -,L  -  1.  It  is  noted  that  to  obtain  the  mapped  value  sk 
corresponding  to  rk,  we  simply  sum  the  histogram  components  from  0  to  rk. 

Local  enhancement  The  histogram  equalization  procedure  is  global  in  the  sense  that 
pixels  are  modified  by  a  transformation  function  which  is  based  on  the  intensity  distribution 
of  an  entire  image.  While  this  global  approach  is  suitable  for  overall  enhancement,  it  is  often 
necessary  to  enhance  details  over  small  areas. 

The  histogram  processing  techniques  are  easily  adaptable  to  local  enhancement.  The 
procedure  is  to  define  an  nxm  neighborhood  and  move  the  center  of  this  area  from  pixel  to 
pixel.  At  each  location  we  compute  the  histogram  of  the  nxm  points  in  the  neighborhood 
and  obtain  a  histogram  equalization  function.  The  function  is  finally  used  to  map  the 
intensity  of  the  pixel  centered  in  the  neighborhood.  The  center  of  the  nxm  region  is  then 
moved  to  an  adjacent  pixel  location  and  the  procedure  is  repeated. 

4.2.3  Thresholding 

Image  thresholding  is  one  of  the  principal  techniques  used  for  object  detection,  especially 
in  applications  requiring  high  data  throughput. 

Suppose  that  the  image  f{x,y)  is  composed  of  a  light  object  on  a  dark  background, 
such  that  object  and  background  pixels  are  grouped  into  two  dominant  modes.  One  obvious 
way  to  extract  the  objects  from  the  background  is  to  select  a  threshold  T  which  separates 
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the  intensity  modes.  Them,  any  point  (x,y)  for  which  f(x,y)  >  T  is  called  an  object  point, 
otherwise  the  object  is  called  a  background.  In  the  case  where  the  image  has  three  dominant 
modes,  we  can  use  the  same  approach  and  classify  a  point  (x,  y)  as  belonging  to  one  object 
class  if  T]  <  /(x,y)  <  T2)  to  the  other  object  class  if  f(x,y)  >  T2,  and  to  the  background  if 
f(x,y)  <  T\.  Based  on  the  foregoing  concepts,  we  may  view  thresholding  as  an  operation 
that  involves  tests  against  a  function  T  of  the  form 

T  =  T[x,y,p(x,y),f(x,y))  (11) 

where  f(x,y)  is  the  intensity  of  point  (x,y),  and  p(x,y)  denotes  some  local  property  of  this 
point,  for  example  the  average  intensity  of  a  neighborhood  centered  at  [x,y).  We  create  a 
thresholded  image  g(x,y)  by  defining 

at*    \      I1  if /(*,»)  >T  ,..,. 

9{J^)=\0  if /(*,y)<T  (12) 

Thus  in  examining  g(x,y),  we  find  that  pixels  labeled  1  correspond  to  objects,  while  pixels 

labeled  0  correspond  to  the  background. 

When  T  depends  only  on  f{x,y),  the  threshold  is  called  global.   If  T  depends  both  on 

f{x,y)  and  p(x,y),  then  the  threshold  is  called  local. 

4.3  Region-Oriented  Segmentation 

The  objective  of  segmentation  is  to  partition  an  image  into  regions.  In  our  case,  we 
are  interested  in  segregating  areas  on  the  hull  of  the  ship  that  need  to  be  cleaned  due  to 
excessive  fouling.  Such  areas  will  be  identified  from  the  images  obtained  using  the  sonars. 

Let  R  represent  the  entire  image  region.  We  may  view  segmentation  as  a  process  that 
partitions  R  into  n  subregions,  /?i,  i?2,  ■■  ■  ,RU,  such  that 

1.  UJLj  Ri  =  R 

2.  #,  is  a  connected  region,  t  =  1, 2,  •  •  • ,  n. 

3.  Rt  n  Rj ■  =  <p  for  all  i  and  j,  i  ^  j 

4.  P{Ri)  =  TRUE  for  »  =  1,2, --.n 

5.  P(Ri  nRj)  =  FALSE  for  i  ?  j 

where  P{Rt)  is  a  logical  predicate  defined  over  the  points  in  set  R^,  and  $  is  the  null  set. 

Condition  1  requires  that  the  segmentation  must  be  complete;  that  is,  every  pixel  must 
be  in  a  region.  The  second  condition  requires  that  points  in  a  region  must  be  connected. 
Condition  3  indicates  that  regions  must  be  disjoint.  Condition  4  deals  with  the  properties 
that  must  be  satisfied  by  the  pixels  in  a  segmented  region.  Finally,  condition  5  indicates 
that  regions  Rz  and  Rj  are  different  in  the  sense  of  predicate  P. 
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Region  growing  by  pixel  aggregation  We  will  start  with  a  set  of  "'seed"  points  and 
from  these  grow  regions  by  appending  to  each  seed  point  those  neighboring  pixels  that  have 
similar  properties  like  intensity.  The  property  that  we  will  use  to  include  a  pixel  in  a  region 
will  be  based  on  the  absolute  difference  between  the  intensity  of  the  seed  and  a  threshold 
value  T. 

Region  splitting  and  merging  The  procedure  discussed  above  grows  regions  starting 
from  a  given  set  of  seed  points.  An  alternative  is  to  subdivide  an  image  into  a  set  of  arbitrary, 
disjoint  regions  and  then  merge  or  split  the  regions  in  an  attempt  to  satisfy  the  conditions 
stated  at  the  beginning  of  the  section.  A  split  and  merge  algorithm  may  be  explained  as 
follows 

Let  R  represent  the  entire  image,  and  select  a  predicate  P.  Assuming  a  square  image,  one 
approach  for  segmenting  R  is  to  successively  subdivide  it  into  smaller  and  smaller  quadrant 
regions  such  that,  for  any  region  R,,  P{R)i)  =  TRUE.  The  procedure  starts  with  an  entire 
region  R.  If  P(R)  is  FALSE,  we  divide  the  image  into  quadrants.  If  P(R)  is  FALSE  for  any 
quadrant  we  divide  that  quadrant  into  subquadrants  and  so  on.  This  particular  splitting 
technique  has  a  convinient  representation  in  the  form  of  a  so-called  "quadtree"  (i.e.  a  tree 
in  which  each  node  has  exactly  four  descendants).  A  simple  illustration  is  shown  in  Fig.3. 
It  is  noted  that  the  root  of  the  tree  corresponds  to  the  entire  image  and  that  each  node 
corresponds  to  a  subdivision. 

If  we  used  only  splitting,  it  is  likely  that  the  final  partition  would  contain  adjacent 
regions  with  identical  properties.  This  may  be  remedied  by  allowing  merging,  as  well  as 
splitting.  In  order  to  satisfy  the  segmentation  conditions  stated  earlier,  we  merge  only 
adjacent  regions  Rt  and  Rj  only  if  =  TRUE. 

The  preceding  discussion  may  be  summarized  by  the  following  procedure  in  which,  at 
any  step,  we 

1.  Split  into  four  disjoint  quadrants  any  region  R*  for  which  P{Rt)  -  FALSE. 

2.  Merge  any  adjacent  regions  R3  and  Rk  for  which  P(Rj  u  Rk)  =  TRUE. 

3.  Stop  when  no  further  merging  or  splitting  is  possible. 

4.4  Description  of  the  sonar  image 

Description  deals  with  the  computation  of  features  like  size,  shape,  etc.,  suitable  for 
differentiating  one  object  from  another  in  the  image.  In  our  case,  we  are  interested  in 
identifying  clusters  of  fouled  areas  and  identify  the  size  and  shape  of  these  clusters. 

Some  simple  descriptors       A  number  of  existing  imaging  systems  are  based  on  regional 
descriptors  which  are  rather  simple  in  nature.  The  area  of  a  region  is  defined  as  the  number 
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of  pixels  contained  within  its  boundary.  The  major  and  minor  axes  of  a  region  are  useful 
for  establishing  the  orientation  of  an  object.  The  ratio  of  the  lengths  of  these  axes,  called 
the  eccentricity,  is  also  a  global  descriptor  of  the  shape.  The  perimeter  of  a  region  is  the 
length  of  its  boundary.  Although  the  perimeter  is  sometimes  used  as  a  descriptor,  its  most 
frequent  application  is  in  establishing  a  measure  of  compactness  of  a  region,  defined  as 

perimeter2  /area. 

Image  description  While  there  are  many  ways  to  describe  an  image,  we  will  primarily 
be  interested  in  the  intensity  distribution  over  a  connected  region  in  the  image.  The  simple 
regional  descriptors  will  provide  us  with  the  shape  and  size  of  the  connected  region.  Subse- 
quently, we  will  use  the  first  and  the  second  moments  of  the  intensities  of  every  pixel  in  the 
connected  region  from  some  fixed  reference  point  to  obtain  the  mean  and  the  standard  devi- 
ation. These  quantities  will  characterize  the  intensity  distribution  over  the  entire  connected 
region  and  will  be  used  to  differentiate  one  connected  region  from  another  in  the  extent  of 
fouling. 
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Figure  5.  A  3  X  3  neighborhood  about  a  point  (x,  y)  in  an  image 
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Figure  6.  A  general  3X3  mask  showing  coefficients  and 
corresponding  image  pixel  locations 
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Figure  7.  (a)  Partitioned  Image,  (b)  Corresponding  quadtree 
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5.0    EXPERIMENTS  WITH  THE  ALUMINUM  PLATE 


5. 1       General 

The  dimensions  of  the  test  plate  used  for  this  study  are  shown  in  Figure  8.  Figure 
6  shows  the  geometry  for  the  test  setup  and  the  relationships  between  slant  range  and 
horizontal  range  for  the  testing  configuration.  Analysis  of  the  geometry  of  Figure  9  would 
indicate  that  if  6  is  the  elevation  angle  of  the  sonar  head  and  the  sweep  azimuth  angle  is  y, 
then  the  horizontal  coordinates  of  a  slant  range  vector  r,  in  global  plate  coordinates  [X,Y,Z 
]  would  be 

X  =  rcos0cosy 

Y  =  rsiny 

Z  =  -rsin0cosy 

In  this  transformation  there  is  a  question  concerning  the  value  to  use  for  the  angle 
6.  Because  the  beam  has  a  height  of  24  degrees,  and  knowing  that  0  is  a  nominal  say  20 
degrees  means  that  the  particular  return  could  be  associated  with  20  plus  or  minus  0/2. 
However,  since  we  are  assuming  that  the  return  is  also  associated  with  the  plane  of  the 
plate,  we  are  able  to  compute  0  in  terms  of  the  particular  slant  range  value,  r,  and  the  height 
value,  h.  It  follows  that 


with  an  uncertainty  associated  with  the  length  of  the  bin  size. 

5.2       Experimental  Data  Files  Obtained 

For  the  experimental  work,  three  basic  data  files  were  obtained  by  sweeping  over 
the  plate  with  the  sonar  set  at  different  gains  and  both  low  and  high  resolution  conditions. 
Too  high  a  gain  setting  would  result  in  filling  the  screen  with  white  dots  while  too  low  a 
gain  setting  would  not  give  a  registration  of  the  roughness  to  be  identified.  The  high  /  low 
resolution  settings  essentially  showed  that  high  resolution  would  be  needed.  The  data  files 
were  labeled  as 

scan_pal_llh.d 

scan_pal_l  ll.d 

scan_pal_5h.d 

The  number  1 1  or  5  refers  to  the  sonar  gain  setting  and  the  designation  1  or  h  refers 
to  the  resolution  of  the  data  to  be  either  low  or  high. 

The  data  from  the  files  have  been  analyzed  by  thresholding,  plotting  on  Cartesian 
coordinates,  contour  plotted,  and,  algorithms  have  been  developed  for  automatic 
identification  of  the  mean  range  to  a  high  spot  as  well  as  effective  size  of  a  region  of  high 
returns. 
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19"  high  -20  degrees 


(0,0)  Y 

Sonar  beam  is  +/- 12  degrees,  1  deg.  wide 


Figure  8  Dimensions  of  the  Aluminum  Test  Plate  Used  for  this  Study 

and  the  Locations  of  the  Nuts  (All  Dimensions  in  Inches) 
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X  =  r  cos8  cosy;  Y  =  r  sin\j/ 


Figure  9  Transformations  from  Slant  Range  to  Plate  Plane  Coordinates 


24 


6.0    RESULTS  OF  CARTESIAN  REPRESENTATION 
OF  THRESHOLDED  RETURNS 

6. 1  General 

Since  the  results  are  contained  in  the  values  of  the  intensity,  it  is  most  convenient  to 
develop  the  data  set  in  terms  of  vector  valued  quantities  where  they  are  entries  in  a  row 
column  matrix.  Intensity  values  are  thus  viewed  as  the  entries  in  the  matrix  Z  (i,  j): 

I=Z(i,j); 

The  index  i,  corresponds  to  the  bearing  where  1  represents  the  initial  heading  of  the 
sonar  and  the  highest  value  represents  the  bearing  as  scan  number  increases.  There  is 
clearly  a  mapping  between  i  and  the  heading  in  degrees.  Similarly,  there  is  a  mapping 
between  j  and  the  range,  j  is  the  bin  number  and,  in  this  case,  takes  on  values  from  0  to 
64.  The  mapping  between  j  and  the  range  in  cm.  is  given  as  previously  noted. 

With  a  view  from  Figure  3  of  the  geometry  of  the  setup,  we  convert  the  slant  range 
data  into  equivalent  x,  y,  coordinates  for  comparison  with  the  Cartesian  plate  coordinates  of 
simulated  fouling  sections. 

The  result  is  shown  in  Figure  10.  Figure  10  contains  the  locations  in  plate 
coordinates  of  the  nuts  that  were  used  to  simulate  discrete  fouling,  and,  superimposed  are 
the  sonar  returns  in  the  thresholded  bands  T=  12  -  15; 

Since  there  is  a  discretization  on  slant  range  of  +/-  9.375  cm.  and  an  angle 
resolution  of  +/-  1  degrees  arc,  it  is  apparent  that  there  are  indeed  returns  that  would  be 
useful  to  work  with,  even  though  the  returns  do  not  lie  precisely  on  the  locations  of  the 
nuts.  The  overall  result  is  quite  remarkable  however,  and  is  very  encouraging. 

Lower  intensity  returns  are  found  over  a  large  area  as  shown  in  Figure  1 1  which, 
for  example  shows  returns  in  the  range  of  intensities  from  7-11. 

A  contour  plot  of  the  matrix  Z  with  thresholded  returns  between  T=[  12,15],  is 
shown  in  Figure  12.  Here,  again,  it  is  clear  that  the  sonar  is  able  to  detect  areas  of 
simulated  roughness.  The  contour  plot  capability  of  the  MATLAB  software  package  has 
been  used  to  produce  the  Figures  for  these  results. 

6.2  Effects  of  Thresholding 

Thresholding  the  data  may  be  done  to  eliminate  the  confusion  cause  by  massive 
amounts  of  low  intensity  data.  Use  of  an  upper  threshold  is  sometimes  useful  if  areas  of 
high  intensity  need  to  be  leveled  for  the  purpose  of  computing  size  of  the  high  intensity 
zones. 

To  threshold  the  data,  the  following  algorithm  is  used. 

for  i=l,m 
for  j= 1,1 

if  Z(i,j)<Tiower.  Z(i,j)  =  Tl0wer; 
if  Z(i,j)>  Tupper,  Z(i,j)  =  TUpper 
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Figure  10  Returns  in  the  Range  12-15  Superimposed  on  the  Simulated  Roughness 
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Figure  1 1  Returns  in  the  Range  7-11  Intensity 
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Fi^e  12.  Contour  plot  of  the 


intensity  matrix  Z. 
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where  Tupper  and  Tiowerare  the  upper  and  lower  thresholds  used  and  m  and  1  are  the  array 
sizes  for  the  number  of  bearings  and  range  bins. 

6.3       Real  Time  Data  Scan  Line  Processing 

In  view  of  the  above,  it  would  seem  that  processing  the  sonar  returns  after  data 
acquisition  could  be  done,  displaying  the  results  for  the  robot  user,  but,  on  the  other  hand, 
with  a  certain  amount  of  on  board  processing,  an  estimation  of  the  location  and  severity  of 
the  fouling  is  really  needed  in  real  time.  In  other  words,  examination  of  the  matrix  Z  (i,  j), 
reveals  that  using  the  intensity  weighted  centroid  computation  would  yield  the  mean 
distance  sO  and  hence,  using  the  bearing  information,  the  mean  locations  xO,  yO,  in  plate 
coordinates  could  be  computed  automatically. 

6.3.1    Mean  Range  Algorithm 

Under  the  assumption  that  the  intensity  /  range  along  a  scanline  is  unimodal,  we  can 
find  the  mean  range  by  the  first  moment  of  intensity.  That  is,  an  intensity  weighted  mean 
range  will  give  a  measure  to  the  'centroid'  of  the  intensity  versus  range  plot.  This 
computation  is  applied  scan  by  scan  and  is  performed  by  the  following  code. 


for  i=l,m 

r0(i)=r*Z(i,:)/(sum(Z(i,:)); 

end; 


where  r0(i)  represents  the  mean  range  along  scan  line  i  for  a  total  number  of  m  bearings. 

6.3.2    Depth  Estimation  Algorithm 

A  measure  of  the  depth  extent  of  the  segments  of  strong  returns  along  a  scan  line 
can  be  found  by  the  second  moment  of  intensity,  and,  having  established  values  for  rO(i), 
we  find  that,  if  a(i)  is  viewed  as  the  standard  deviation  of  intensity  weighted  variations 
from  the  mean  range,  then  the  algorithm  is  as  follows. 

for  i=l,m 

c2(i)=(r-r0(i)*ones(size(r)).A2*Z(i,:)  /  (sum(Z(i,:)); 

a(i)=sqrt(a2(i)); 
end; 


In  the  above,  .A2  refers  to  the  element  by  element  square  of  the  vector  of  range 
differences  from  the  mean  range  for  that  bearing. 

6.3.3    Mean  Strength  Algorithm 

The  mean  strength  of  returns  contained  within  the  area  of  the  assumed  unimodal 
high  strength  zone  is  also  an  interesting  metric  and  is  found  by  averaging  the  intensity 
values  containd  inside  the  ranges  defined  by  (r0(i)  -  a(i))  and  (r0(i)  +  c(i)).  We  call  this 
metric,  the  mean  strength  (mst(i)). 
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Such  a  compulation  yields  the  result  shown  in  Figure  13.  Notice  that  the  locations 
[xO,  yO],  do  not  precisely  correspond  to  the  locations  of  the  highest  intensity  returns 
because  they  take  into  account  the  total  effect  of  weighting  all  intensities  within  the 
thresholded  data  set. 

Processing  the  data  to  estimate  the  depth  variance,  yields  Figure  14  which  shows 
that  only  the  central  cluster  of  returns  has  significant  linear  extent.  While  this  idea  is  able  to 
estimate  points  of  roughness,  it  will  be  necessary  to  test  the  algorithms  on  data  that  do 
come  from  real  roughness  with  some  significant  extent  along  bearing  lines. 

The  estimation  of  the  mean  strength  of  returns  along  a  scan  line  coming  from 
returns  that  have  significant  extent  is  given  by  the  data  in  Figure  15  where  the  means 
strength  (mst)  is  plotted  versus  the  bearing  angle.  Most  of  the  returns  do  not  exceed  the 
lower  threshold,  and  so  the  mst  is  zero.  However,  those  that  are  associated  with  a 
significant  return,  show  the  mst  to  be  between  12  and  15  as  expected.  A  lower  mst  is  an 
indicator  that  more  than  one  return  is  being  used  in  the  averaging  as  shown  by  Figure  16, 
and  that  the  distribution  is  not  unimodal.  The  average  mst  is  lower  but  is  still  of  interest 
and  could  be  useful. 

What  the  authors  contend  is  that  the  hull  cleaning  Robot  will  eventually  need  to  be 
given  an  [xO,  yO]  coordinate  together  with  a  significant  width  and  depth  measure  and  an 
intensity  measure  such  as  the  mst,  that  could  be  computed  automatically  scan  by  scan  . 
Further  work  to  test  the  concept  is  needed. 
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Figure  13         XO,YO  as  Computed  from  Processed  Returns 
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Figure  14  Depth  Estimation 
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Figure  15  Mean  Strength  Estimation 
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Figure  16  Indicauon  of  Multi  Modality 
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CONCLUSION 


In  conclusion,  an  initial  study  has  begun  to  identify  whether  or  not  sonar  scans  of  ship 
plating  could  be  a  viable  method  of  determining  the  location  and  extend  of  hull  fouling.  A 
simulated  set  of  barnacles  using  one  half  inch  steel  nuts  has  been  ensonified  and  processing 
of  returns  from  an  ST  725  sonar  in  one  configuration  has  shown  that  the  locations  of 
significant  roughness  could  indeed  be  found.  What  is  more,  the  algorithms  presented  are 
expected  to  lead  to  a  real  time  capability  for  an  autonomous  robot  to  locate  and  analyze  the 
sections  of  hull  fouling.  A  necessary  extension  of  this  work  is  to  investigate  the  ability  of 
the  sonar  and  methods  presented  here  to  be  successful  against  real  marine  growth. 
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